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Abstract 

Ring topologies of repressing genes have qualitatively 
different long-term dynamics if the number of genes is 
odd (they oscillate) or even (they exhibit bistability) . 
However, these attract ors may not fully explain the 
observed behavior in transient and stochastic envi- 
ronments such as the cell. We show here that even re- 
pressilators possess quasi-stable, travelling- wave pe- 
riodic solutions that are reachable, long-lived and ro- 
bust to parameter changes. These solutions underlie 
the sustained oscillations observed in even rings in 
the stochastic regime, even if these circuits are ex- 
pected to behave as switches. The existence of such 
solutions can also be exploited for control purposes: 
operation of the system around the quasi-stable orbit 
allows us to turn on and off the oscillations reliably 
and on demand. We illustrate these ideas with a sim- 
ple protocol based on optical interference that can in- 
duce oscillations robustly both in the stochastic and 
deterministic regimes. 

Key words: synthetic biology; oscillatory gene 
networks; generalized repressilator; traveling waves; 
stochastic dynamics 

Introduction 

Recent experimental advances in cellular and molec- 
ular biology have made it possible to engineer intri- 
cate gene regulatory circuits (1). Inspired in many 
cases by electronic elements, simple gene networks 



have been designed to perform reproducible, low-level 
functions. Some classic examples include the toggle 
switch ([2]), the genetic ring oscillator known as the 
repressilator (3), or a circuit that can exhibit both 
oscillatory and switching behavior through the alter- 
ation of biochemical interactions Such simple cir- 
cuits could be potentially interconnected and built- 
up to form more elaborate 'biological devices' with 
large numbers of components. This trend is facili- 
tated by simulation software containing large number 
of genes ([5]) as well as libraries of composable biologi- 
cal parts for experimental realization (|5J). Simple syn- 
thetic modules can also be integrated into the com- 
plex machinery of the cell, as in the oscillator recently 
implemented in a mammalian cell ([7J), or interfaced 
with cellular pathways to induce particular responses, 
as in the construct where the toggle switch was con- 
nected to the SOS pathway to induce DNA protection 
mechanisms in E. coli when exposed to UV light (|SJ). 
Similar principles have been exploited in the rational 
design of internal negative feedback operated in con- 
junction with external arabinose-driven positive feed- 
back to produce cell-synchronized oscillations (j9]). 

The central role played by oscillations in cellu- 
lar function has made oscillatory circuits a primary 
target for the analysis and design of synthetic net- 
works. A particular area of interest is the elucidation 
of strategies leading to robust timing and sequen- 
tial activation in the cell. For instance, key stages 
in developmental biology and in cell differentiation 
may be controlled by so-called master regulators — a 
small set of transcription factors sequentially acti- 
vating and driving several other genes with accurate 
timing (fTUl [TT | [T2 ]) . addition, studies of both natu- 
ral ([T3l fl4j) and engineered circuits ([T5]) indicate that 
the correct timing and order of gene activation is a 
key characteristic of balanced, optimal cell function, 
as it reduces the metabolic burdening that ensues 
from the continuous presence of heterologous sub- 
stances (fT6|). 

In this paper, we consider the dynamics and con- 
trol of noisy genetic oscillatory circuits in quasi-stable 
mode operation. We exemplify our results with one of 
the simplest and most widely studied synthetic net- 
works: the TV-gene ring repressilator (Fig. la). Some 
natural networks of master regulators ([TO]) contain 
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such ring structures as subnetworks, making the ex- 
ploration of their dynamic behavior relevant for both 
naturally occurring and synthetic systems. The un- 
derlying idea is well-known: when observing the dy- 
namics of systems operating in highly variable envi- 
ronments, such as the cell, it might not be enough 
to characterize only the long-term attractors of the 
system since unstable solutions can play a significant 
role. For instance, quasi-stable transients might be 
so long-lived as to be the most significant feature of 
the observed noisy dynamics (fT7j) . Moreover, the 
presence of noise in nonlinear systems may induce 
non-stationary dynamics in systems with only fixed 
point attractors in the deterministic setting ([T8j) or, 
conversely, noise may act as a stabilizer of unstable 
deterministic states ([T9]) . 

In the generalized repressilator, results due to 
Smith ([201 |2T|) show that rings with an even num- 
ber of genes (e.g., the toggle switch (2) with N = 2) 
exhibit multistability and hence behave like switches 
in the stochastic regime. On the other hand, rings 
with an odd number of genes (like the standard re- 
pressilator (j3j with N = 3) have a globally attract- 
ing limit cycle and are therefore oscillators both in 
the stochastic and deterministic regimes. However, 
we show here that generalized repressilators possess 
an intricate structure of unstable periodic orbits that 
play an important part in their observable noisy and 
transient dynamics. In particular, even rings have a 
quasi-stable limit cycle which, although unstable in 
terms of linear Floquet stability analysis, has only 
one unstable direction with a very slow escape rate. 
This means that trajectories are attracted to the limit 
cycle from all directions but one, hence leading to 
long-lived, inducible periodic transients in the deter- 
ministic setting and to sustained oscillations in the 
stochastic system. These effects become more pro- 
nounced as the number of genes grows. Therefore the 
finite-time, observable noisy dynamics of an even re- 
pressilator ring is not necessarily static (switch-like) 
but rather exhibits oscillatory characteristics. 

In addition to their effect on the observable dy- 
namics, quasi- stable oscillatory modes can be used as 
operating points to control the system around them. 
The advantage of such a scheme is that the oscil- 
lations can be switched on and off, unlike limit cy- 



cle attracting behavior. Operation around unstable 
modes, usually illustrated with the example of the 
inverted pendulum ([22]) . is a classic scenario in con- 
trol theory for enhanced controllability and speed of 
response. It has a long and successful history of ap- 
plications in fluid flow control ([23)) and in the steering 
of jet aircraft ([24]) . Here, we illustrate the applica- 
tion of this concept to gene networks with a simple 
protocol of controlled interference based on an opti- 
cal mechanism for readout and induction of gene ex- 
pression. Our simulations show that even repressila- 
tor rings in quasi-stable operation behave as a robust 
and on-demand switchable oscillator in which genes 
become upregulated periodically in an ordered se- 
quence according to a travelling-wave solution. This 
characteristic, which is robust at high and low copy 
numbers, could be used for synthetic biological ap- 
plications such as accurately timed interference with 
naturally occurring networks. 

Theory 

Model equations and stability of fixed points 

The generalized repressilator consists of a ring of N 
genes in which transcription of each gene is repressed 
by the product of the preceding gene (Fig. la). A 
deterministic model of this circuit is given by the fol- 
lowing set of ordinary differential equations: 

m = — — 2 c 2m 

Pi = c 3 mi-C4Pi, (1) 

where pi and rrii describe protein and mRNA con- 
centrations for each gene, respectively (|3j) . Here, 
i = 1,...,7V with the periodic boundary condition 
Po — Pn, and c\ (cs) is the creation rate and C2 (04) is 
the degradation rate for the mRNAs (proteins). The 
production of mRNA is modeled as a source term 
that depends nonlinearly on the concentration of the 
inhibitor protein. Proteins are assumed to be pro- 
duced at a rate linearly dependent on the amount of 
the corresponding mRNA. The degradation of mRNA 
and proteins is assumed to be linearly proportional 
to their current amount. The toggle switch ([2]) and 
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the repressilator oscillator (3), which have both been 
implemented in E. coli, are special cases with N = 2 
and N = 3, respectively. Based on analytical results 
on monotone systems due to Smith ([2Q| |2TL the sta- 
bility analysis of this circuit reveals a fundamental 
difference between rings with odd and even numbers 
of genes. We briefly sketch some of the main results, 
which are also summarized in Fig. lb. 

The stability analysis characterizes the long-term 
dynamic behaviors of the deterministic system. An 
example of such behavior is given by the fixed points 
of the system, i.e., the states in which the dynam- 
ics is stationary. The variation of a parameter can 
produce a change in the stability or the existence of 
fixed points or other attractors. This is called a bi- 
furcation and it leads to qualitative changes in the 
long-term behavior of the system. One can find the 
parameter values at which bifurcations are produced 
by performing a bifurcation analysis, which can be 
carried out analytically (in some simple cases) or nu- 
merically with the aid of continuation software pack- 
ages such as AUTO ([25]). which can also track the 
stability of periodic solutions. 

In our system ([!]), the fixed points, where all 
derivatives are zero, are found from the condition: 

p*(l+pU 2 ) = — = c, Vi. (2) 

The parameter c defined in Eq. ([2| will play the role 
of the bifurcation parameter for even rings. A posi- 
tive and real solution is obtained if all proteins have 
the same concentration: p* = = p m , Mi and 
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This solution, which exists as long as c is positive, 
is stable for small c and becomes unstable for larger 
values of c in both odd and even rings. 

In the case of even rings, a pitchfork bifurcation 
takes place at c — 2 for all N. The two additional 
stable fixed points arising at that value of the pa- 
rameter correspond to p* = p* +2 ^ Vi, which 



gives: 



= \ ~ = Pd- ( 3 ) 

Note that p u — > c — 1/c and pd — > 1/c for large c. 
The new fixed points of the system ([!]) correspond 
to two distinct dimerized states: one in which genes 
with odd indices are upregulated (p u ) while genes 
with even indices are downregulated {pa)\ and an- 
other symmetric state where the genes with odd and 
even indices exchange their patterns of regulation. 
These solutions are equivalent to tiling the ring with 
copies of the steady state solution of the two-gene 
ring, and are similar to other dimerized degenerate 
solutions in classic models of conjugated polymers 
and spin chains ([26]) . Therefore, after the bifurcation, 
the system is bistable, i.e, it behaves like a switch in 
the presence of noise. 

In the case of odd rings, p m becomes unstable fol- 
lowing a bifurcation that occurs at a value c(N) that 
approaches 2 as TV grows. However, in this case the 
bifurcation is Hopf: no additional fixed points appear 
but rather the bifurcation signals the emergence of a 
periodic solution. Smith ([20]) proved that in mono- 
tone systems (i.e., systems in which partial deriva- 
tives do not change sign) such as the repressilator, the 
periodic solution that emerges is a globally attract- 
ing stable limit cycle. Therefore odd rings behave as 
stable oscillators following the Hopf bifurcation. 

Floquet theory and unstable periodic orbits 

The stability analysis presented above does not pro- 
vide information about unstable periodic solutions. 
Although, in principle, unstable periodic orbits are 
not relevant for the long-term deterministic dynam- 
ics, they can be key to the observed dynamics, es- 
pecially if the orbits involve slow time scales. Such 
quasi-stable oscillations can appear as transients in 
deterministic simulations and are likely to be ob- 
served in the corresponding stochastic simulations. 
In fact, it was in numerical simulations that we first 
noticed the relevance of these modes in even repres- 
silator rings. 
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Floquet theory can be used to find periodic so- 
lutions and quantify their linear stability in terms 
of their Poincare map, i.e., the crossings of the or- 
bit with a (hyper)plane in phase space. Under this 
analysis, a periodic solution (a closed orbit) becomes 
a fixed point of the Poincare map and its stability 
is reformulated as the linear stability of this fixed 
point. The eigenvalues of the Poincare map linearized 
around the fixed point constitute the Floquet multi- 
pliers. They indicate how an infinitesimal pertur- 
bation around the orbit decays or grows (exponen- 
tially). The periodic solution is linearly stable if all 
the Floquet multipliers have magnitudes smaller than 
unity (see the Supplementary Material for references 
and details on Floquet Theory). In some cases, a 
few (possibly only one) Floquet multipliers will be 
slightly larger than one. We will then have 'quasi- 
stable' periodic solutions in that it takes a long time 
to diverge away from them. Quasi-stability in this 
sense is a local property. To assess if these solutions 
will be reachable (and therefore observable in the dy- 
namics), one needs to employ global techniques, e.g. 
sampling the space of initial conditions. However, 
Floquet analysis provides an indication of the possi- 
bility of observable, yet unstable, periodic solutions. 
If a periodic orbit has a small number of very weakly 
unstable directions, it is likely that it could be ob- 
served as long-lived periodic transients in the deter- 
ministic system and that it could also play a role in 
the stochastic dynamics. Moreover, such quasi-stable 
orbits are good candidates for a control mechanism 
that can induce oscillations on demand, as shown be- 
low. 



Methods 

Numerical simulations and analysis of the dy- 
namics 

The deterministic system of ODEs ([I]) was solved nu- 
merically with an adaptive fourth-order Runge-Kutta 
integrator ([27|) . in which the step-size automatically 
adapts to meet the required accuracy e. We have 
checked that the inducibility, reachability and tran- 
sient times of the quasi-stable oscillations are not af- 



fected by the accuracy of the integrator by using the 
Runge-Kutta integrator with accuracies e between 
10 -2 and 10 -8 (see Supplementary Material for de- 
tails). In addition, the observability of the unstable 
orbits was confirmed by using a nonlinear integra- 
tor (|28j) . 

The bifurcation analysis and the calculation of the 
Floquet multipliers of the unstable periodic orbits 
were carried out with the numerical continuation soft- 
ware AUTO (25) (see Supplementary Material). 

Stochastic simulations of the generalized repressi- 
lator were performed using the classical Gillespie al- 
gorithm (29). Random numbers and quasi-random 
numbers for numerical simulations were generated 
with the GSL Scientific Library (|30j) . 

Global robustness analysis and control aspects 
of quasi-stable oscillations 

As part of our numerical evaluation of the general- 
ized repressilator, we have developed a method to 
carry out a robustness and reachability analysis of 
its quasi-stable oscillations. This was necessary be- 
cause available global robustness tools (28) quantify 
changes in fixed points induced by parameter vari- 
ations and are therefore not directly applicable to 
oscillations. In order to evaluate the global robust- 
ness and inducibility of the quasi-stable oscillations, 
we attempt to induce sustained oscillations with a 
predetermined intervention and quantify changes in 
the observed response when the model parameters are 
varied. The method defines an operating point in pa- 
rameter space (the reference set c*), based on biologi- 
cally appropriate estimates, and a hyper cube around 
it to account for biological variability, temperature 
gradients and other noise. We then sample param- 
eter sets from the hypercube using Reverse Halton 
Sequences (|3T]h quasi-random sequences that have 
been shown to converge faster than standard Monte 
Carlo sampling for high dimensional spaces (2). For 
each sampled parameter set, we attempt to induce 
oscillations with the STOP-KICK scenario described 
below and record if the system evolves towards sus- 
tained oscillations. If oscillations are observed, we 
calculate numerically the period of the oscillation and 
the change in shape (see supplementary material for 
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details). Characterizing the change in shape is es- 
sential to establish that the oscillation remains de- 
tectable and functionally recognizable in the biolog- 
ical system. Note that here we are only concerned 
with global robustness of the reachability of the so- 
lution. A modification of the same algorithm could 
be used to study the parameter combinations that 
contribute most strongly to the sensitivity of the net- 
work, a question relevant for the experimental tuning 
of the system that is not addressed here. 

Results 

Stable and quasi-stable oscillations in the gen- 
eralized repressilator 

As pointed out in the Theory section, odd repressila- 
tor rings are globally attracted to stable limit cycle 
oscillations for c > 2. Numerical simulations show 
that the period of these solutions increases linearly 
with the number of genes in the ring (Fig. 2d). The 
stability analysis also shows that, in contrast, even 
rings only support fixed points as stable solutions. 
However, direct dynamical simulations of even repres- 
silator rings reveal the existence of long-lived periodic 
solutions, which are easily reachable, as checked by 
extensive sampling (not shown) of the space of initial 
conditions. The period of these oscillations also in- 
creases linearly with the number of genes, albeit with 
a slope that is approximately half of that in odd rings 
(Fig. 2a). 

These numerical observations do not pose a con- 
tradiction with the stability analysis above: the ob- 
served oscillations in even rings are periodic solu- 
tions yet unstable. Unstable solutions can be stud- 
ied using the numerical bifurcation detection software 
AUTO ([25]) . a continuation package that does not 
rely on dynamical simulations (see Supplementary 
material). We have used AUTO to find bifurcations 
in the biologically relevant range of the parameter c 
and to assess the linear stability of fixed points and 
periodic solutions — the latter through Floquet anal- 
ysis. 

The result for even rings is presented in Table 1. 
In agreement with the analytical results, a pitchfork 



bifurcation is found numerically as a branching point 
at c = 2, above which Hopf bifurcations leading to 
the appearance of unstable periodic solutions are de- 
tected in all rings with more than 4 genes. The Flo- 
quet stability analysis indicates that the first unstable 
periodic orbit to emerge has only one unstable direc- 
tion, regardless of the number of genes. The only 
positive Floquet multiplier, which indicates how fast 
the trajectory diverges away from the orbit, is small 
and decreases as the length of the ring increases. This 
is the signature of quasi- stability: if this periodic or- 
bit is reached, it will be long-lived. We have also 
checked that this solution is reachable through nu- 
merical sampling of the space of initial conditions 
(see Supplementary material). Such reachable quasi- 
stable modes affect significantly the observed tran- 
sient dynamics and also play a central role in stochas- 
tic dynamics where unstable solutions are explored 
under the effect of noise. Both these conditions are 
relevant for dynamics of genetic circuits inside the 
cell. 

The existence of quasi-stable modes provides us 
with the opportunity to design a distinct control 
strategy. If we control the system to revolve around 
a quasi-stable mode, the result is an oscillator that 
can be switched on, kept oscillating and switched off 
on demand. Below, we introduce a simple implemen- 
tation of such a scenario and evaluate its robustness 
of operation. Note that an intricate family of unsta- 
ble periodic orbits with high symmetry exists both 
in odd and even rings (see Table 1 and Supplemen- 
tary Material). However, these additional periodic 
solutions have several unstable directions that make 
them essentially unobservable and uncontrollable. 

Spatio-temporal structure of the periodic so- 
lutions 

The spatio-temporal structure of the periodic solu- 
tions, both in the odd and even cases, corresponds 
to a travelling- wave solution propagating around the 
ring. The snapshots in Figs. 2b, c show that this prop- 
agation occurs against the backdrop of the dimerized 
fixed point solution of the even ring, where all odd 
(even) numbered genes are 'up' while the even (odd) 
numbered genes are 'down' pi). Clearly, a dimerized 
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configuration cannot be accommodated in an odd 
ring. This leads to a kink-like (frustrated) solution 
where two consecutive genes have similar expression 
levels. This local imbalance of repression induces a 
dynamical instability that makes the kink propagate 
around the odd ring in a periodic fashion (Fig. 2b). 
This spatio-temporal structure underlies the limit cy- 
cle solution in odd rings. The fact that the period of 
the limit cycle increases (roughly) linearly with the 
number of genes indicates that the speed of propaga- 
tion of the kink is (roughly) constant. 

The quasi-stable periodic solution in even rings can 
be interpreted under the same prism. Figure 2b shows 
that it corresponds to two interacting kinks propagat- 
ing around the ring at a roughly constant speed with 
a period that is approximately one half of that of the 
closest odd ring (Fig. 2a). The instability of this pe- 
riodic solution has a clear meaning in this picture: if 
the two kinks 'collide', they annihilate each other and 
the system returns to the stable fixed point, i.e., the 
dimerized solution. Figure 2b also shows that each 
kink has a minimum spatial width that depends on 
the parameters of the model. Hence it is more dif- 
ficult to find these oscillatory solutions in rings that 
are not large enough to fit two such perturbations al- 
though they can still be observed in smaller, biolog- 
ically realizable rings (see Supplementary material). 
For clarity, we have chosen to illustrate the spatio- 
temporal structure of the solutions with long rings. 
However, we have checked that the quasi-stable pe- 
riodic orbits in rings with N = 6, 8, 10 (not shown) 
maintain the features of the two-kink-structure and 
operate under the same principles as the long rings 
shown in Fig. 2b. 

The spatio-temporal structure of the periodic solu- 
tions in repressilator rings shows a strong parallelism 
with similar dynamical solutions observed in classical 
models of discrete lattices (33). The travelling- wave 
nature of the oscillations could have potential biolog- 
ical applicability if one were to use this circuit as a 
control element for genes that must be activated in a 
particular order and for a pre-defined time interval. 



Robust induction of quasi-stable oscillations in 
the deterministic regime 

The causal constraints imposed by the travelling- 
wave structure means that the manipulation of one 
gene will have a predictable effect on the others. This 
property can be used to induce and stop oscillations 
in even rings reliably by activation of one gene for a 
short time span. We illustrated this simple scenario 
in Fig. 3a. First, the even ring is forced to converge 
to one of the fixed point solutions with a STOP sig- 
nal that consists of the external activation of gene j 
for a time interval longer than the period of the oscil- 
lation. This signal is used to 'initialize' the system, 
suppressing any transient oscillations present in the 
system. Once the system is at rest, the oscillation 
can be started with a KICK signal, consisting of the 
external activation of gene j + 1 with a step function 
of width and amplitude similar to those of the oscil- 
latory pattern. Such signals can be imparted non- 
invasively through an optical mechanism that uses 
UV or red light to activate the production of mR- 
N As of particular genes (j34[ [35j) . 

In order to check that the proposed protocol is ro- 
bust to parameter variations, we have carried out a 
global robustness analysis as outlined in the Meth- 
ods section. We construct a hypercube by taking 
variations of 5%, 10% and 20% around the refer- 
ence values of the parameters in Eq. and take 
10 4 samples in this hypercube varying all parameters 
simultaneously. Sampling is performed with quasi- 
random Reverse Halton Sequences for improved con- 
vergence (2, 31). Figure 3b shows that the fraction of 
parameter samples that lead to oscillations with this 
protocol converges to 1 for large rings. Our numer- 
ics also show that oscillations can be elicited with 
significant robustness in rings with N > 6. When 
oscillations are present, the period shows very small 
variation with respect to the reference set, as shown 
by the coefficient of variation (Fig. 36, inset). We 
have also quantified the change in the shape of oscil- 
lations through a normalized mean square measure 
and found that the shapes in the perturbed system 
exhibit very high similarity (« 99%) to the reference 
set (see Supplementary material). In summary, our 
global robustness analysis indicates that in the deter- 
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ministic regime, long-lived quasi-stable periodic solu- 
tions are reliably inducible in larger rings with a sin- 
gle intervention. For moderate size rings, oscillations 
can still be induced for a large fraction of the param- 
eter hyper cube but are shorter lived. This suggests 
that repeated interventions could be used in order to 
keep the ring in the quasi- stable oscillating state. A 
simple control protocol that implements these ideas 
is proposed in the following section and shown to be 
applicable for rings as small as N — 6 operating in 
the stochastic regime. 

Stochastic oscillations in even rings and 
readout-based control 

We have used the standard Gillespie algorithm ([29|) 
to study the generalized repressilator in the stochas- 
tic regime, i.e., when intrinsic noise is high due to 
low copy numbers. It is well known that stochastic 
models of odd rings behave as oscillators and that the 
travelling wave structure is preserved (3, 36). In the 
case of even rings, we have performed long stochas- 
tic simulations (not shown) and found bistability and 
switching events, as expected from the long-term at- 
tractors of the underlying deterministic system. Ad- 
ditionally, the simulations show sustained oscillatory 
behavior, especially in longer rings although they are 
also observable in rings as small as N = 6 (see Sup- 
plementary material). The oscillations appear in a 
variety of settings: as transients from a variety of ini- 
tial conditions; spontaneously emerging from one of 
the stable points; or associated with switching events. 
It is also easy to induce such oscillations with local- 
ized interventions in particular genes; to sustain them 
with periodic driving; and to terminate them with a 
prolonged induction of a gene (as in the STOP sig- 
nal above). We have examined the structure of these 
oscillations and they correspond well with the quasi- 
stable periodic solutions in the deterministic system: 
their period is approximately half the period of the 
closest odd ring (Fig. 2 a) and the spatio-temporal 
travelling- wave structure is maintained. 

Our numerical simulations confirm the relevance 
of the underlying quasi-stable oscillations for the ob- 
served stochastic dynamics of even rings. Similarly 
to the deterministic case, the quasi-stable mode can 



also be used as a control operating point such that the 
system becomes switchable. The robust reachability 
of this mode allows us to use an extremely simpli- 
fied feedback mechanism that could be implemented 
through an optical read-out (GFP, YFP or luciferase 
protein labeling) and response (on-demand UV or 
red light gene transcription activation) ([34]). The 
simple control scheme illustrated in Fig. 4 a uses the 
optical read-out from two successive proteins in the 
ring to introduce optical KICK signals that sustain 
the oscillation based on a threshold rule (Fig. Ab). 
The oscillation can be started and terminated using 
the same optical signals. Although we have chosen 
to illustrate the possible implementation of the con- 
trol scheme with light sensitive inducers, it is worth 
remarking that any suitable mechanism capable of 
precisely-timed gene induction with good spatial res- 
olution at the cell population could be used. A poten- 
tial advantage of this switchable mode of operation 
is the economical and targeted use of the transcrip- 
tional resources without overburdening the cell with 
unnecessary mRNA production ([T6]) . 

Discussion 

In this work, we have studied how the presence of 
quasi-stable periodic solutions affects the observable 
dynamics of even repressilator rings. Previously, even 
rings have been thought of as switches due to the fact 
that they only support fixed point solutions. How- 
ever, our bifurcation analysis reveals the existence 
of a set of unstable orbits, some of which have slow 
timescales associated with them. These quasi-stable 
periodic solutions are both reachable and long-lived, 
thus playing a role in the observed dynamics, both 
transient and stochastic. This suggests that oscil- 
latory behavior might be more widespread than ex- 
pected in genetic models since it could feature in sys- 
tems that possess only static attractors. 

The presence of quasi-stable solutions provides us 
with the possibility of designing control protocols 
that operate the system around such modes so that 
the oscillations can be turned on and off reliably. 
Our numerics indicate that a robust mechanism could 
be implemented based on appropriate optical feed- 
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back to switch the system between stable fixed points 
and quasi-stable oscillations. Although the proposed 
pared-down control scheme is only intended to pro- 
vide an illustration of the potential implementation 
and its performance could be improved with the use 
of optimized strategies for stochastic and robust con- 
trol that take into account specific details of the ex- 
perimental setup, it is worth discussing some of its 
limitations. The challenge for the dynamical control 
algorithm is to deliver the optical interference sig- 
nal necessary for the induction of gene expression for 
a short time period, to a particular spatial area of 
the cell population, and with a well-controlled delay 
following the fluorescent signal of the proteins. The 
proposed scheme shows both enough spatial and time 
resolution to address individual cells in a popula- 
tion with well-defined pulses (|34|) . The scheme would 
need to rely on proper calibration of life-times of flu- 
orescent proteins affected by phototoxic and photo- 
bleaching effects, as reviewed in detail by Bennett 
and Hasty ([37]). Finally, since the delay between the 
actual protein concentration and the corresponding 
fluorescent signal introduced by the maturation time 
(~ 2-8 min) is short compared to the period of the 
oscillation (~ 130 min), this would be acceptable for 
the control scheme. 

A synthetic circuit operating under such principles 
could be interfaced with a naturally occurring net- 
work to induce an intrinsic interference that is inter- 
ruptible on-demand. The switchability of this regula- 
tory element can help avoid the appearance of adverse 
cumulative effects. The NF/^B pathway is an example 
where such a regulator could provide controlled acti- 
vation over short time intervals as an alternative to 
conventional knock-downs and other functional inter- 
ventions that modify the balance of important pro- 
teins for the cell cycle (38 , 39). The underlying travel- 
ling wave structure of the observed periodic solutions 
could also be potentially useful for design purposes. 
It allows for coordinated intervention when the tim- 
ing and order of activation of different pathways is 
crucial. Examples of cellular networks, e.g. in de- 
velopmental biology, indicate that timed patterns of 
sequential activation are at the heart of the function 
of families of master regulators (fTT| [14| \42 \ |43 |) and, in 
the case of the vertebrate segmentation clock ([4Q| [4T ]) . 



the associated oscillations are well-defined but do not 
survive in the long-term. The importance of hetero- 
geneously timed gene induction has also been high- 
lighted in a model of Arabinose uptake in E. coli ([T3]). 
Experiments with genetically engineered yeast have 
also shown that pulsed activation of chaperons fol- 
lowed by pulsed activation of the associated heterol- 
ogous proteins is more efficient at maximizing the 
production of particular metabolites (fT6|) . These ap- 
plications hint at potential uses for circuits that can 
produce sequential patterns of activation on demand, 
such as the even repressilator studied here, which in- 
teract with other cellular pathways via intrinsic pro- 
teins, thus avoiding the timed delivery of external 
agents through the cell membrane. 

The design of control strategies for the operation 
of systems around an inherently unstable state has a 
long history in other disciplines (e.g., flight and fluid 
control) since it affords enhanced responsiveness and 
controllability with relative ease and simplicity of de- 
sign (22). This strategy differs fundamentally from 
the biochemical alteration of the network topology 
proposed by Atkinson et al. (j4j) based on a smaller 
gene circuit but with a complex regulatory scheme 
involving promoter and repressor sites regulating one 
gene. The molecular kinetics of such regulators are 
less well understood than those with single regulatory 
sites due to unavoidable cross-talk and compound 
logic. The ring topology relies on simple regulation 
to provide a sequence of causal signals but at the 
expense of involving a larger number of genes. 

The present scheme is also in contrast with previ- 
ously engineered gene circuits, such as odd repressi- 
lators, which possess globally attracting limit cycles 
leading to behavior that is robust yet not control- 
lable. Quasi-stable operation, on the other hand, is 
robustly switchable. The switchability of the oscilla- 
tor coupled with dynamic control that affords good 
spatial resolution could be used to elicit localized os- 
cillations in cell populations as an aid to examine 
mechanisms of cell synchronization. It is an open area 
of current research to elucidate the role of a design 
concept based on control around unstable behavior, 
similar to the inverted pendulum in classic control 
theory, to further our understanding of cell strate- 
gies and its potential use in the design of synthetic 
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topologies that can interfere with naturally occurring 
pathways. 
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Bifurcation analysis and unstable periodic solutions of repressilator rings with even number of genes. 

We use the continuation package AUTO (|25[) to obtain the bifurcations of rings of size N . The parameter c, defined 
in Eq. |2|, is swept in the biologically relevant range c G [0.001,30] by changing c\ with C2 = 0.12, C3 = 0.16 and 
C4 = 0.06 constant. In agreement with analytical calculations, a branching point (P), corresponding to a pitchfork 
bifurcation, is found at c = 2. A series of Hopf bifurcations (HB) linked to the emergence of unstable periodic 
solutions are found subsequently. Floquet analysis indicates that the first unstable orbit to emerge has only one 
unstable direction, regardless of the dimension of the system, and that the maximal Floquet multiplier decreases 
with increasing N. Hence this periodic solution is quasi-stable: if it is reached, the divergence away from it is slow, 
and gets slower for longer rings. Other unstable orbits are present but their high instability makes them irrelevant 
to the observed dynamics. A similar structure of unstable orbits exists in odd rings (see Supplementary Material). 
The figure on the right shows the bifurcation diagrams for even rings of length N = 6, 12, 16. The unstable periodic 
orbits, shown in red dashed lines, emerge through Hopf bifurcations. 



Figure Legends 

Figure [TJ Attractors of the generalized repres- 
silator model 

(a) Topology of the generalized repressilator: N genes in a 
cycle where each gene is repressed by the protein product 
of the preceding gene. Also shown is the reaction scheme 
underlying the dynamical system |l]) with production and 
degradation terms for the mRNA (rrij) and protein (pj) 
of each gene. The repression of the production of mRNA 
is modelled by a Hill- type term H(pj-i). (b) Typical 
time traces of the long-term deterministic dynamics of 
an odd ring and an even ring above the bifurcation point 
c = 2: odd rings converge to a globally attracting periodic 
solution while even rings converge to fixed points. The 
time traces shown correspond to N = 23 and N = 22. (c) 
Stability of the fixed points of the system as a function 
of the bifurcation parameter c. Even rings undergo a 
pitchfork bifurcation at c = 2, leading to the emergence 
of two stable fixed points. Odd rings undergo a Hopf 
bifurcation leading to the emergence of a limit cycle. The 
critical parameter for the Hopf bifurcation depends on 
N but tends to c = 2 as N grows (see Supplementary 
material) . 

Figure [2} Periodic solutions and travelling 
waves in the generalized repressilator model 

(a) The period of the limit cycle of the deterministic 
model of odd rings (solid line) increases linearly with the 
length of the ring. Simulations of the stochastic version 
of this system using the Gillespie algorithm show that the 
period (shown as circles) follows the same trend, although 
they are slightly larger. The period of the quasi-stable so- 
lutions found in even rings (deterministic and stochastic) 
increases linearly with the number of genes but with a 
slope that is half that of the odd rings. The inset shows 
representative time traces of the periodic solutions in odd 
rings (stable) and in even rings (quasi-stable), (b) Time 
snapshot of the spatial distribution of the concentrations 
of two successive proteins concentrations for the periodic 
solution in the odd ring with N = 23. The solution has 
a travelling-wave structure with a kink-like perturbation 
propagating around the ring, indicated by the arrow in 
the bottom figure. The bottom figure represents the min- 
imum distance |Apj| = min(|p w — Pj\, \pd — Pj\) between 



the traveling wave solution and the dimerized solution 
with an alternating pattern of protein expression given 
by p u and pa- The distance becomes large around the 
kink in the traveling wave solution, (c) Same as (b) for 
the quasi-stable periodic solution of the even ring with 
N — 22. In this case, the travelling wave solution has two 
kinks that propagate around the ring, as indicated by the 
arrows. 

Figure [3} Robust induction of oscillations in 
even rings in the deterministic regime 

(a) The quasi-stable periodic solution in even repressila- 
tor rings can be induced with a simple sequence of sig- 
nals. First, apply a STOP signal to gene j to force the 
system to approach a fixed point solution. Second, ap- 
ply a KICK signal to gene j + 1 to drive the ring into 
oscillation. The oscillation can be terminated at will by 
applying another STOP signal. The signals can be im- 
plemented via on-demand UV or red light gene transcrip- 
tion activation ((34}. This STOP-KICK-STOP protocol 
is shown here for a ring with N = 18 with parameters 
Cl = 2.6, c 2 = 0.12, c 3 = 0.2, c 4 = 0.06. (b) Global 
robustness analysis of the inducibility of the quasi-stable 
oscillations. The STOP-KICK scenario is applied to 10 4 
random combinations of parameters c» for each even ring 
of length N and we record the proportion of parameter 
sets that lead to five oscillations. The parameters are 
sampled with reverse Halton sequences from a hypercube 
with 5%(B), 10%(«) and 20%(4) variation around the ref- 
erence set. Quasi-stationary oscillations are robustly in- 
duced for N > 10, while smaller rings can be kept in 
the oscillating state applying repeated interventions in a 
simple control protocol as shown in Fig. [4] (Inset) The 
oscillations are robust in shape (not shown) and in pe- 
riod to changes in the parameters. The relative variabil- 
ity (coefficient of variation) of the period of the induced 
oscillations is small and decreases with the length of the 
ring. 

Figure [4[ Stochastic oscillations in even rings 
and readout-based control. 

(a) Illustration of the readout-based control scheme for 
a ring of 6 genes. Two proteins of the ring are read out 
with fluorescent tags. This readout is then compared to a 



reference denned according to the oscillating behavior of 
the ring, with similar period and a shift between consec- 
utive genes. The reference comparison is threshold-based 
and leads to an ON-OFF (1-0) control for the KICK sig- 
nals. These can be implemented with light responsive 
genes promoters. In the numerical simulations shown in 
(b), the KICK signals are indicated with the red markings 
in the upper panels, (b) A simple readout-based control 
reliably switches on the oscillations, sustains them and 
switches them off. The control mechanism functions by 
monitoring two successive proteins in the ring. When- 
ever each of them falls below a threshold, a KICK signal 
for the corresponding protein is given. These threshold- 
based KICK signals are indicated with red and magenta 
markings in the upper panels. The oscillation can be 
terminated with a STOP signal as in the deterministic 
state. The optical read-out can be based on GFP or YFP 
protein labeling while the response can be implemented 
with on-demand UV or red light that enhances the pro- 
duction of the corresponding mRNAs (|34p . The figure 
shows the application of this mechanism to a ring with 
N — 10. The stochastic time traces correspond to the 
protein expression of proteins pj with j = 1, 3, 5, 7, 9 and 
the corresponding control (top) in response to proteins p\ 
and p2 ( trace not shown). The right figure is a magni- 
fication of the dashed square inside the main figure. We 
have also checked that this control protocol is applicable 
for rings with as low as N = 6 genes (not shown). 
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1 Stability analysis and bifurcations in odd and even rings 
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Table SMI: Bifurcation analysis of repressilator 
bifurcation diagrams. 



5 3 Gene Ring 




with odd number of genes and corresponding 




Figure 5: Bifurcation diagrams for odd and even case. The unstable periodic orbits are shown in red dashed 
lines. As the number of genes in the ring increases the first unstable orbit moves towards stable structures 
(stable fixed points for the even rings and stable limit cycles for the odd rings. This suggests that the first 
unstable orbits may become more observable in larger rings. The sampling of initial condition space shown 
below confirms this hint.) 



Table SMI and Figure [5] summarize the result of our analysis using AUTO to obtain the bifurcations 
of odd rings of size N. As in the main text, the parameter c is swept in the biologically relevant range 
c G [0.001,30]. In agreement with analytical calculations summarized in the main text, a series of Hopf 



bifurcations (HB) are found. The first one leads to a stable limit cycle, whose maximum Floquet multiplier 
is 1 corresponding to the trivial direction along the orbit. The subsequent bifurcations give birth to unstable 
periodic orbits with an increasing number of unstable directions (always even). These unstable orbits are 
essentially irrelevant to the observed dynamics because they are highly unstable and not reachable. The 
relationships between the periods of oscillations in rings of different lengths hints at a strong symmetry in 
the spatio-temporal structure of these solutions. We summarize the findings in Table SMI and Figure [5j 

2 Basic Floquet Theory and Quasi-Stable Oscillating Orbits under 
Noise 




Figure 6: Poincare map with quasi-stable fixed point (black) and possible trajectories. The periodic orbit 
(red trajectory) is linearly stable along the left diagonal, where the arrows are pointing towards the fixed 
point and unstable along the right diagonal, where the arrows are pointing away from the fixed point. The 
escape scenario from unstable periodic orbit due to pertubations is illustrated with black trajectory. 

As pointed out in the main text we have observed long lasting oscillations in the even rings randomly 
sampling the initial condition space. The observed oscillating modes can be assessed regarding their stability 
with the Floquet Multipliers theory ([3]). Here, the question about stability of the orbit can be reformulated to 
the question for the stability of the corresponding fixed point x* in the Poincare map. If vq is an infinitesimal 
perturbation such that x* + is in 5, then after the first return to S 

x* + v 1 = P(x* + v ) « P(x*) + [DP(x*)]v (4) 

The Eigenvalues Xj of the linearized Poincare map [DP(x*)] obtained this way are nontrivial Floquet mul- 
tipliers. The closed orbit is linearly stable if and only if |Aj| < 1 for all j. We illustrate the case where the 
closed orbit has unstable directions and how the system may escape from it in Figure [6] above. The Poincare 
section (shown in gray) has stable directions (left diagonal) and unstable directions (right diagonal). The 
system may escape the closed orbit (black trajectory) in a directed manner once perturbed along the unstable 
direction. 

Since it is a biological system it is important to show that the quasi-stable orbits are observable under noise 
effects. We can simulate small biological noise effects using less accurate integration with larger maximal 



error e per step (see "Methods" section in the main text). This is in addition to the robustness analysis 
under parameter variations illustrated in Figure SSI. It turns out that the oscillations do not become less 
robust and last as long as with the accurate integration. 




10 15 20 



Figure SSI: a Robustness of STOP-KICK-STOP scenario as illustrated in the main text Figure 3a. The 
accuracy control parameter e for circles is 10 -2 , whereas in the main text the parameter was 10 -8 . Both 
curves are for the same total variation. Comparing the result to the main text, we see that less accurate 
integration does not change the robustness of the STOP-KICK-STOP scenario significantly. Less accurate 
integration is a way to simulate the biological noise and these quasi-stable oscillations do not become signifi- 
cantly less robust under the noise simulated with less accurate integration as shown here and also simulated 
with Gillespie algorithm as shown later in the control scheme in the main text. (See the paragraph below 
for the intuitive explanation for this phenomenon.) b Different accuracy also does not significantly affect 
the length of transient oscillations, which die off within a certain period of time (for instance within 50 os- 
cillations). The figure illustrates the mean lengths of oscillations, which occur if starting from the randomly 
chosen initial condition and inspect the region for 50 oscillations. The error bars are first and third quartiles. 



We try to give an intuitive explanation for the observed effect. Errors resulting from less accurate inte- 
gration affect all variables describing protein and mRNA concentrations. It has the effect that the system is 
perturbed in all directions due to this noise. However, the generalized repressilator system can escape the 
oscillating orbit through directed movement along one out of many possible directions (see Figure [6] above 
and Table 1 in the main text). If the system is perturbed along other directions, it is pulled back to the 
limit cycle, because all other directions are attactive. Therefore random perturbations do not neccessarily 
kill the transient oscillations. On the contrary, random perturbations may interfer with the directed escape 
out of the limit cycle and this might even stabilize the oscillations. In fact in the following section we show 
fluctuations due to internal noise stabilize the oscillations. 



3 Sampling of Initial Condition Space and Spontanious Transient 
Oscillations in the Stochastic Setting 



The bifurcation analysis revealed the existence of quasi-stable periodic solutions in the generalized repressi- 
lator model. However, it is not clear if these solutions are observable under random sampling of the initial 
condition space. Here we show for the deterministic setting how often on average a randomly sampled initial 
condition can lead to at least 5,10, and as many as 50 oscillations Fig. [7] 
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Figure 7: Sampling of the initial condition space for the deterministic simulations. The initial condition 
hyper cube defined as to maximal amplitude (Note: the maximal amplitude of the proteins is not the 
same as the maximal amplitude for the mRNAs) has been sampled 10 4 times using pseudo random numbers 
(reverse Halton sequences). 



In the stochastic setting the noise due to small copy numbers does not destroy the oscillations in the 
even rings. On the contrary, this is also how we have discovered these modes. The quasi-stable structure 
identified with the deterministic bifurcation theory becomes even more prevalent under internal noise. So, in 
this case noise does actually excite the oscillations. The following parameters have been used c\ — 1.6, = 
0.12, cs — 2.6, C4 = 0.06. The picture shows that in the stochastic regime the system spontaneously goes into 
the oscillating state more frequently than the equivalent deterministic description. The initial conditions for 
the stochastic simulation are chosen to be the same as in the deterministic case and 10 runs for each starting 
value has been collected. 



Typical transient oscillating state in 6-Gene -Ring 
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Figure 8: Transient oscilations are more prevalent in the stochastic setting. For the same initial conditions 
and the same constants the system goes more frequently into the oscillating state than in the deterministic 
setting. 



4 Robustness Analysis and Shape Similarity 

We perform global robustness analysis, which means all parameters are varied at the same time within 
specified multidimensional hypercube. If we want to obtain a global robustness of let say 10 gene ring, 
then the dimension of the according parameter space would be 40, there are 4 constants for each gene 
(ci, C2, C3, C4) and these are varied independently of each other. This procedure is fundamentally different to 
standard robustness assessments, where only one parameter is gradually varied at the same time. 

The parameter hypercube is created as 5%, 10%, 20% variation around the reference state, where all 
parameters are the same for all the genes. We sample this parameter hypercube in an efficient manner 
using pseudo-random numbers created with reverse Halton sequences. This sampling of large dimensional 
space with pseudo-random numbers was shown to be more efficient than the classical MC (j2j). After selecting 
a set of parameters in this way, we evolve the system in time first applying the STOP signal, which puts 
it into stationary up/down solution. Then we apply a KICK scenario and record if the oscillations have 
persisted for at least 5 periods. The requirement here is that all genes have shown a non decaying amplitude. 
In this way we obtain 10 4 samples for each gene number and plot the fraction of the samples leading to the 
oscillations on the y-Axes. 

In addition to the 0,1 decision on if or not the fraction oscillates we quantify the change in the shape 
of the oscillation. This is independent of the variation in period and amplitude. The measure ssh would 
give large dissimilarity scores if the oscillations would become spiky, the gene being in the up state is much 
smaller than being in the down state or vise versa. It would give intermediate scores if instead of being 
rather quadratic, it becomes like a Gaussian shape. Mathematically, we first scale the time for the reference 
shape and for the perturbed system shape to 1. Through the scaling we get rid of the variations in the 
period lengths. Then we apply the standard normalized Z^-norm, 

SSh =< 7T re f(t)7T pert (t) > I < 7T ref (t) > / < 7T pert (t) > (5) 

where the normalization scales away the differences in the amplitude. This measure assesses shape changes 
for moderate parameter perturbations and gives an additional information to the variation in the oscillation 
amplitude and period length. 
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Figure 9: Shape variations under the global parameter changes. On the left we see the numerically observed 
oscillations, which differ moderately in the period length and amplitude. The reference shape for the unper- 
turbed system is shown in red. On the right the time scale (x-Axis) has been scaled with the period and the 
maximal amplitude for each curve respectively. This demonstrates how the scaling in time and normalized 
L2 score assesses the shape similarity. 



5 Bifurcation Software AUTO 

AUTO (|4]) performs systematic variation of specified parameters and assessment of qualitative characteristics 
of ODE system such as fixed points, periodic orbits and their stability. The ODEs are of the form: 

u'(t) = f(u(t),p), f(;-),u(-)eR n (6) 

The limited bifurcation analysis is done by solving algebraic equations 

f(u,p)=0, f(-,.),ueR n (7) 

One of the limitations of the package is that the continuation must start nearby guessed/known stable fixed 
point. If such starting point is given, when solution branches are calculated as one or more parameters are 
gradually changed. 

For our system there is a regime, where only one stable fixed point exists (see the Theory section in the 
main text). Starting from that equilibrium point we have performed the continuation gradually changing 
one parameter and identified the bifurcation points. 
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